COMPUTING THE FIRST EIGENPAIR FOR 
PROBLEMS WITH VARIABLE EXPONENTS 

MARCELLO BELLOMI, MARCO CALIARI, AND MARCO SQUASSINA 

Abstract. We compute the first eigenpair for variable exponent eigenvalue problems. We compare the 
homogeneous definition of first eigenvalue with previous nonhomogeneous notions in the literature. We 
C^ highlight the symmetry breaking phenomena. 

o 

^ 1. Introduction 

■^^ The aim of this paper is the study of numerical sohitions to the minimization problem introduced in [7] 

O n l^ ^ _ =„. l|V?/||p(.) 



o 



(1.1) Ai = inf 



•^ Here il is a bounded domain in M" and the variable exponent p : Cl ^ IR+ is a smooth function such that 

1 < p~ < p{x) < p"*" < oo for every a; G f2. The norm || • Wptx) is the so-called Luxemburg norm 



(1-2) ||/|lp(,)-inf{7>0 



p{x) I 

p{x) ~ . 



G If p is a constant function, the problem reduces (up to a power p) to the minimization of the quotient 

^ (1.3) inf 

(^ and, as it is known, its associated Euler-Lagrange equation is 
00 



n 



a 



(1.4) -div(|Vu|P^2VM)=A|M|^"^M. 

In this case, we refer the reader to [12, 13] for the theoretical aspects and to [2] for a recent numerical 
^r^ analysis. The special case p = 2 is the classical eigenvalue problem for the Laplacian — Au = Au, for 

T-H which we refer the reader to [11]. In general, in these type of problems, is crucial that some homogeneity 

►^ holds, namely, if w is a minimizer, so is uju for any non-zero real constant lo. On the contrary, the quotient 



H (1.5) J^ 

^ / |zi|P(-) 

Jo. 

with variable exponents fails to possess this feature. Therefore, as we point out in the following section, 

its infimum over nontrivial functions of Wf^ ^ {Vt) turns out to be often equal zero and no minimizer 

exists [4,6]. A way to avoid this collapse is to impose the constraint ^^\uY''^^^dx = C. Unfortunately, 

doing so, minimizers obtained for different normalization constants C are difficult to compare. For a 

suitable p(x), it could even happen that any A > is an eigenvalue for some choice of C. Thus (1.5) is 

not a proper generalization of (1-3), which has well defined (variational) eigenvalues, although the full 

spectrum is not completely understood yet. A way to avoid this situation is to use the Rayleigh quotient 
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(1.1), restoring the necessary homogeneity. In the integrand of (1-2), the use of the measure p{x)~^dx 
simplifies the Euler-Lagrange equation. The Sobolev inequahty [3] ||w||p(a;) < C'llVw||p(2,), shows that 
Ai > 0. It is easy to see that (1.1) has a non-negative minimizer. Pick a minimizing sequence of Vj, 
namely ||fj||p(a;) = 1 and |lV'i;j||p(a;) = Ai +o(l). By Relhch theorem for variable Sobolev exponents [3], 
up to a subsequence, we find u such that Vj -^ u in LP*^^^(il) and Vvj -^ Vu in LP*^^^(ri). This yields 



Ai < ||Vu| 



p{x)/ 



\u\\p{x) < limj^ool 



Vw 



'j\\pix)/\\V3\\p(x) 



— Ai. Notice that if u is a minimizer so is \u\ > 0. 



By the maximum principle of [9] , u has a fixed sign. In [7] the Euler-Lagrange equation for a minimizer 
is derived. Precisely, it holds 



(1.6) 

where we have set 

K = ||Vu| 

More generally, A € 

(1.7) 



Vu 



K 



p(x). 



p{x)-2 



(^,Vr;)=Ai5 



p{x)-2 



k 



V, V77 e c^{n), 



k = 



\p{x). 



Ai = 



K 



is eigenvalue if there exists u € Wq (^), i* ^ 0, such that 



S = 



l,p(x). 



p(x) 



Vu 



K 



p(x) 



Ik 



p(x) 






V77) = AS* / 
Jo. 



p{x)-2 u 



r;, ^ri^C^in). 



It follows from the regularity theory developed in [1] that the solutions to (1.7) are continuous provided 
that p{x) is Holder continuous. If Ai is the minimum in (1.1), we have A > Ai in (1.7), thus Ai is the first 
eigenvalue and a corresponding solution is the first eigenfunction. Contrary to the constant exponent 
case [12,13], it is currently unknown if, in the variable exponent case, the first eigenvalue is simple, and 
if a given positive eigenfunction is automatically a first one. Concerning higher eigenvalues, in [15] the 
authors have recently proved that there is a sequence of eigenvalues of (1.7) with \j /" 00 and if 

1 1 ■ 

then there are constants Ci, C2 > 0, that depend only on n and p^, such that 

(1.8) Ci |f7| A"/(i+'") < #{j : A^ < A} < C2 \^\ A"/(i-'") for A > large, 

where |fi| is the Lebesgue measure of fi. Observe that, in the case of constant p, (1.7) reduces not exactly 
to (1.4), which is homogeneous of degree p — 1, but rather to the problem (homogeneous of degree 0) 



<1, 



iVujP-^Vi 



= A 



ip-i 



u^Wl'^iSl) 



^ iiv^iir' 

Thus (1.8) should be compared to Ci \^\ A"/^ < #{j : A^ < A} < C2 \^\ A"/^, obtained in [8]. 

1.1. A different notion in the literature. We compare the minimization procedure with the Rayleigh 
quotient with Luxemburg norm and that without it, namely 



|Vw|P(^) 



a 



inf 

In this framework, if A <E R and u S Wq (^) then (u, A) is called eigenpair if m 7^ and 

|Vu|P(")-2^Vu, Vr/) = A /" |u|P(")^2y?7, V77 € I^o^^^"^(r!). 

Set A = {A > : A is an eigenvalue}. It is well known [12, 13] that, if the function p(a;) is constant, then 
the problem has a sequence of eigenvalues, sup A = +00 and inf A > 0. In the general case, it follows 



/ 
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from [6] that A is a nonempty infinite set and sup A ~ +00. Define A* := inf A. We recall that we often 
have A* = (recall that Ai > in (f.f)). Consider the following Rayleigh quotients 



inf 






M* 



inf 



IVwp(^) 






l.P(-),r« / ly|p(x) 



Then, in [6], the authors prove that A* > <=> /^* > <;=^ /I^ > 0. Furthermore, if there is an open subset 
[/ C ri and a point a;o € C/ such that /'(xq) < p(x) (or >) for all x € dU , then A* = [6, Theorem 3.1]. 
In particular, if p(a;) has strictly local minimum (or maximum) points in fi, then A* = 0. There are also 
statements giving some sufficient conditions for inf A > 0. Let n > 1. If there is a vector i g M" \ {0} 
such that, for any a; e fi, the map t M- p{x + t£) is monotone on {t : x + t£ ^ ft}, then A* > [6, Theorem 
3.3]. If n == 1, then A* > if and only if the function p{x) is monotone [6, Theorem 3.2]. 




Figure 1. Unit square, p{x,y) = 5 + 3sin(37ra;), a;-axis (left), y-axis (right) and top views 



2. An algorithm to compute the first eigenpair 

In this section we briefly describe an algorithm to approximate Ai in (1.1) and compute the corre- 
sponding eigenfunction. We start defining 



Ai 



inf 



|V^||^(.) 



R{u) 
uew^-"'-'''' {a)\{o} \W\\l(x) uew^-'"'''\n)\{o} S{u) 



inf 



A? 
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Figure 2. Unit disk, p{x,y) = 11 + 9sin(27ra;), a;-axis (left), y-axis (right) and top views 



It is now possible to apply the inverse power method, where the (j + l)-th iteration is 
(2.1a) 



vP^ — argniin(i?(w) — VS'(u-')u) ~ argmin J(u) 



(2.1b) 
(2.1c) 



„j+i - 



A^' 



5(uJ + l)l/2 



where u^ is the result of the previous iteration and, by (2.1b), has Luxemburg norm equal to 1. It is 
possible to show (see [5]) that the algorithm converges to a critical point of R{u)/S{u), even if it is not 
possible in general to prove convergence to the smallest eingevalue. However, a good choice of the initial 
guess M° can reasonably assure that the result is the smallest eigenvalue. For the computation of \7S{u)'ri, 
for given u and r], in the so called inner problem (2.1a), we observe that ii u ^ its Luxemburg norm 
7(u) — ||u||p(a.) is implicitly defined by 



(2.2) 



F(u,7) 



u{x) 



1 



p{x) 



p{x) 



1 = 



Therefore, we can use the differentiation of implicit functions to get 

p{x)-2 

VuF[u,-i)T] Jn 



VII 



\\p(x)V 



u 

\u\\p{x) 



\p{x) 



d^F{u,-f) 



\\p{x) 



p{x) 
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Figure 3. Annulus, p(x,y) — 4 + 2sin(27rx), x-axis (left), j/-axis (right) and top views 



from which 



WS{u)tj = 2 



Mpix) 



p{x)-2 



wq 



WW 



p{x) 



p(x) 



Since we are mainly interested is some particular two-dimensional domains, such as a rectangle, a disk 
or an annulus, we approximated the problem by the finite element method which well adapts to dif- 
ferent geometries by constructing an appropriate discretization mesh. In particular, we used the tool 
FreeFEM-l— I- [10] which can handle minimization problems as (2.1a) through the function NLCS (nonlinear 
conjugate gradient method, Fletcher-Reeves implementation). Such a function requires the application 
of the gradient of J(u) to a test function 77 

p(x)-2 I. _, p(x)-2 

(Vm, Vr;) 



yj{u)r] 



Vu 

\^^p(x) 



h^lpix) 



w'rj 



Vu 



iVul 



p{x) 



P{x) 



hmpix) 



P{x) 



and an initial guess which, for the (j + l)-th iteration of the inverse power method, is Uj/A\. The stopping 
criterion for the inverse power method is based on the difference of two successive approximations of Ai . 

2.1. Some details of the algorithm. The algorithm described above requires recurrent computations of 
the Luxemburg norm of a function. For a given m 7^ 0, it is the zero of the function F{u, •) defined in (2.2). 
This is a C^(0, +00) convex and monotonically decreasing function in 7, with lim^_j.o+ F(u,j) = +00 



and lim, 



■7— v+00 



F(m,7) 



-1. Therefore, it is possible to apply the quadratically convergent Newton's 



method in order to find its unique zero, starting with an initial guess 70 on its left hand side (i.e., such 
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that i^(u,7o) > 0). As pointed out above, the inverse power method cannot guarantee the convergence 
to the smallest eigenvalue and relative eigenfunction. It is very reasonable to expect that if the initial 
guess u" for the method is close enough to the smallest eigenfunction, then the algorithm will converge 
to it. For p = 2 the problem essentially reduces to the Helmholtz equation 

—Au = Xiu, 

for which the eigenfunctions are well known for the domains we have in mind. Therefore, starting with 
p = 2 and u'^ the eigenfunction corresponding to the smallest eigenvalue for the Helmholtz equation we 
moved to the desired p{x) through a standard continuation technique. 

We tested our algorithm using both linear and quadratic finite elements and checked the correct order 
of convergence (two and three for the rectangle and one and two for the circular domains, respectively). 
Moreover, we checked that the results with constant p were consistent with those reported in [2]. The 
results in the next section were obtained with quadratic finite elements. We observed convergence of 
our algorithm also for some cases with p{x) < 2, for which the Hessian of J{u) degenerates and the 
nonlinear conjugate gradient is not guaranteed to converge. In this case some authors add a regularization 
parameter to the functional J{u) (see [2]). We moreover observed sometimes slow convergence of the 
nonlinear conjugate gradient. In this case, a more sophisticated method, using the Hessian of J{u) or 
an approximation of it, could be employed. Another possibility would be to use a preconditioner, either 
based on a low order approximation of J{u) (see again [2] and reference therein) or on a linearized version 
of J(u). The implementation of a more robust and fast algorithm, on which we are currently working, is 
beyond the scopes of this paper. 




Figure 4. p{x,y) = 5 + 3sin(37ra;), ^(a;, y) = 11 + 9 sin(27ra;) &nAp{x,y) = 4 + 2sin(27r2;), respectively 



2.2. Examples and breaking symmetry. We show in this section the results we obtained by applying 
the algorithm described above to three test cases, in the square, the disk and the annulus, respectively. 
For each test, we report the plot of the obtained eigenfunction from the x-axis, the y-axis and from the 
top view, respectively. 

The first case (Figure 1) refers to the unit square [0, 1] x [0, 1] and p{x, y) = 5 + 3 sin(37ra;). The first plot 
(left) in this figure does not present special features because the exponent p{x,y) is independent of y. In 
the second plot (right) the profile is different as it feels the diffusion variation in the x-variable. Both 
plots are symmetric with respect to the center of the domain, since p{x,y) has a center of symmetry in 
(1/2, 1/2), see Figure 4. The second case (Figure 2) refers to the unit disk with center and radius 1 
and p{x,y) = 11 + 9sin(27ra;). The maximum of p{x,y) is quite high (see Figure 4) and the profile is 
reminiscent of the one for the limiting case p = oo. The second plot is remarkable because we loose the 
symmetry for the center of the domain. The third case (Figure 3) refers to the annulus with center 0, 
external radius 1 and internal radius 0.25 and p{x,y) —4 + 2sin(27r2:). The resulting eigenfunction is 
more shifted than the case in the unit square and the case in the unit disc. Even the shape of the domain 
influences the contour, here we see that the annulus reflects intensely the mold of the exponent. The first 
plot (left) maintains the center of symmetry in (0,0) and p{x,y) does not depend on y. 
Already in the one dimensional case, it is evident that the the logarithm of first eigenfunction is not a 
concave function, in general, contrary to the constant exponent case where this was proved to be true [16]. 
See Figure 5 for an example in this regard. 
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Figure 5. Interval [—1, 1], p{x) = 28 + 26cos(27ra;), eigenfunction (left) and its logarithm (right) 

In the case where n has some symmetry and p{x) is a radially symmetric (resp. axially symmetric with 
respect to some half-space) function, then some symmetry (resp. partial symmetry) results where recently 
obtained in [14] for semi-stable solutions and mountain-pass solutions. 
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